function K = bending_bb(dofs,V,T,d,desc)
%        K = bending_bb(dofs,V,T,d,desc)
% This function computes the bending matrix $\int_\Omega \Delta \phi_i \Delta \phi_j$. 
% and this is build in quick version for convenient access.
% here is the reference element
%                V2
%                |  \
%                |    \
%                |      \             
%                |        \
%                |          \
%                |            \
%                |              \
%                |                \
%                V3---------------- V1

% prepare the global bending matrix for saving time
nt = size(T,1);
% get the areas of all triangles.
x21 = V(T(:,2),1) - V(T(:,1),1);
x31 = V(T(:,3),1) - V(T(:,1),1);
y21 = V(T(:,2),2) - V(T(:,1),2);
y31 = V(T(:,3),2) - V(T(:,1),2);
areas = (x21.*y31 - x31.*y21)/2;

Mat = build_dim2(d-2);  Id = diag(ones((d+1)*(d+2)/2,1));
Du = d*de_cast_step(Id,d,1,0,-1,desc); % the direction derivate on v1-v3
Duu = (d-1)*de_cast_step(Du,d-1,1,0,-1,desc); % the direction derivate on v1-v3
Duv = (d-1)*de_cast_step(Du,d-1,0,1,-1,desc);  % blending partial derivatives
Dv = d*de_cast_step(Id,d,0,1,-1,desc);  % the direction derivate on v2-v3
Dvv = (d-1)*de_cast_step(Dv,d-1,0,1,-1,desc); % the direction derivate on v2-v3
Int_UU = Duu'*Mat*Duu; Int_UV = Duu'*Mat*Dvv;
Int_UW = Duu'*Mat*Duv; Int_VV = Dvv'*Mat*Dvv;
Int_VW = Dvv'*Mat*Duv; Int_WW = Duv'*Mat*Duv;
        
m = (d+1)*(d+2)/2; m2 = m*m*nt;
Indx1 = zeros(m2,1); Indx2 = zeros(m2,1);
S = zeros(m2,1); pos = 1;
for k = 1:nt
    V1=V(T(k,1),:);V2=V(T(k,2),:);V3=V(T(k,3),:);
    u = (V2(2)-V3(2))^2+(V2(1)-V3(1))^2;
    v = (V1(2)-V3(2))^2+(V1(1)-V3(1))^2;
    w = 2*((V2(2)-V3(2))*(V1(2)-V3(2))+(V2(1)-V3(1))*(V1(1)-V3(1)));
    LocK = (u*u*Int_UU - u*w*(Int_UW+Int_UW') + v*v*Int_VV +...
        w*w*Int_WW - v*w*(Int_VW+Int_VW') + u*v*(Int_UV+Int_UV'))/(16*areas(k)^3);   % J = 2*tri_area(K)
    [i,j,s] = find(LocK);
    L = length(i);
    Indx1(pos:(pos + L-1)) = dofs(k,i)';
    Indx2(pos:(pos + L-1)) = dofs(k,j)';
    S(pos:(pos + L-1)) = s;
    pos = pos + L;   
end
dim = max(max(dofs))-min(min(dofs))+1; pos = pos -1;
K = sparse(Indx1(1:pos),Indx2(1:pos),S(1:pos),dim,dim);